Lyapunov exponents in constrained and unconstrained ordinary differential equations 
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We discuss several numerical methods for calculating Lyapunov exponents (a quantitative measure 
of chaos) in systems of ordinary differential equations. We pay particular attention to constrained 
systems, and we introduce a variety of techniques to address the complications introduced by con- 
straints. For all cases considered, we develop both deviation vector methods, which follow the 
time-evolution of the difference between two nearby trajectories, and Jacobian methods, which use 
the Jacobian matrix to determine the true local behavior of the system. We also assess the merits 
of the various methods, and discuss assorted subtleties and potential sources of error. 

PACS numbers: 05.45.Pq, 05.45.-a, 95.10.Fh 



I. INTRODUCTION 

Chaos exists in a wide variety of nonlinear mathe- 
matical and physical systems, and ordinary differential 
equations are no exception. Since the original discov- 
ery by Edward Lorenz of deterministic chaos in a toy 
atmosphere model (consisting of twelve differential equa- 
tions) 0, a seemingly endless variety of ODEs exhibit- 
ing extreme sensitivity to initial conditions has emerged. 
Many tools, both qualitative and quantitative, have been 
developed to investigate this chaotic behavior. Perhaps 
the most important quantitative measure of chaos is the 
method of Lyapunov exponents, which indicate the av- 
erage rate of separation for nearby trajectories. (See 
0- C3- d 0- E3- for some recent investigations into mea- 
sures of chaos and their applications.) The present paper 
is concerned with general methods for calculating these 
exponents in arbitrary systems of ODEs. We first re- 
view the techniques for calculating Lyapunov exponents 
in unconstrained systems 0, (where each coordinate 
represents a true degree of freedom) , and then introduce 
several new methods for calculating Lyapunov exponents 
in constrained systems (where there are more coordinates 
than there are degrees of freedom). 

A defining characteristic of a chaotic dynamical sys- 
tem is sensitive dependence on initial conditions, and the 
Lyapunov exponents are a way of quantifying this sensi- 
tivity. In a system of ordinary differential equations, this 
sensitive dependence corresponds to an exponential sep- 
aration of nearby phase-space trajectories: if two initial 
conditions are initially separated by a distance eo, the 
total separation grows (on average) according to 



e(t) = e e 



A/. 



(1.1) 



where A is a positive constant (with units of inverse time) 
called the Lyapunov exponent. Two important caveats 
to Eq. (|l.lfl are necessary. First, this prescription yields 
only the largest Lyapunov exponent, but a dynamical 



'Electronic address: mhartl@tapir.caltech.edu 



system with n degrees of freedom has in general n such 
exponents. Second, Eq. Ql.lf) does not constitute a rigor- 
ous definition, since it defines a true Lyapunov exponent 
only if e is "infinitesimal." A more precise definition of 
Lyapunov exponents involves the true local behavior of 
the dynamical system, i.e., the derivative or its higher- 
dimensional generalization. 

We can go beyond Eq. l|l.lfl to determine (at least in 
principle) all n Lyapunov exponents by considering not 
just one nearby initial condition, but rather a ball of ini- 
tial conditions with radius eo- As discussed in Sec. [n] 
this ball evolves into an n-dimensional ellipsoid under 
the time-evolution of the flow, and the lengths of this 
ellipsoid's principal axes determine the Lyapunov expo- 
nents. We will see that there are many advantages to 
this ellipsoid view, both conceptual and computational. 

We discuss in Sees. |H] and IIIII several techniques for 
calculating Lyapunov exponents in ODEs, and compare 
the relative merits of the various methods. We take 
special care to explain methods for the calculation of 
all n Lyapunov exponents. Our principal examples are 
two well-studied and simple systems: the Lorenz equa- 
tions (Sec. HID II) and the forced damped pendulum 
(Sec. Ill D 2|) . The techniques and code were developed 
and tested on the much more complex problem of spin- 
ning bodies orbiting rotating (Kerr) black holes, as dis- 
cussed briefly in Sec. IIII Dl and at length in |10L fll| - 

Our two model systems are unconstrained, so that 
each variable represents a true degree of freedom. As 
we see in Sec. IIIII following the evolution of a phase- 
space ellipsoid — and hence calculating the Lyapunov 
exponents — becomes problematic when the system is 
constrained. Such systems are common in physics, with 
constraints arising for both mathematical and physical 
reasons. For example, instead of using the angle 9 to de- 
scribe the position of a pendulum, we may find it math- 
ematically convenient to integrate the equations of mo- 
tion in Cartesian coordinates [x, y), with a constraint 
on the value of x 2 + y 2 . Another example is a spin- 
ning astronomical body, whose spin is typically described 
by the components of its spin vector S = (S x , S y , S z ). 
On physical grounds, we might wish to fix the magni- 
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tude ||S|| = S = JS% + SI + 5|, so that only two of the 

three spin components represent true degrees of freedom. 

We describe in Sec. Mil three methods for finding Lya- 
punov exponents in constrained systems. Our principal 
example of a constrained system is the forced damped 
pendulum described in Cartesian coordinates, a system 
chosen both for its conceptual simplicity and to facilitate 
comparison with the same system without constraints. 
We also show the application of these techniques to the 
dynamics of spinning compact objects in general relativ- 
ity. It was the investigation of these constrained systems 
in that led to the development of the key ideas de- 
scribed in this paper. 

We have developed a general-purpose implementation 
of the principal algorithms in this paper in C++, which 
is available for download 01 ■ The user must specify the 
system of equations (and a Jacobian matrix if necessary), 
as well as a few other parameters, but the main proce- 
dures are not tied to any particular system. Most of the 
results in this paper were calculated using this implemen- 
tation. 

We use boldface to indicate Euclidean vectors, and 
the symbol log signifies the natural logarithm log e in 
all cases. We refer to the principal semiaxes of an n- 
dimensional ellipsoid as "axes" or "principal axes" for 
brevity. 

II. LYAPUNOV EXPONENTS IN 
UNCONSTRAINED FLOWS 

There are two primary approaches to calculating Lya- 
punov exponents in systems of ordinary differential equa- 
tions. The first method involves the integration of two 
trajectories initially separated by a small deviation vec- 
tor; we obtain a measure of the divergence rate by keep- 
ing track of the length of this deviation vector. We re- 
fer to this as the deviation vector method. The second 
method uses a rigorous linearization of the equations of 
motion (the Jacobian matrix) in order to capture the true 
local behavior of the dynamical system. We call this the 
Jacobian method. Though computationally slower, the 
Jacobian method is more rigorous, and also opens the 
possibility of calculating more than just the principal ex- 
ponent. In this section we discuss these two methods, 
and several variations on each theme, in the context of 
unconstrained dynamical systems. 

When discussing Lyapunov exponents in ordinary dif- 
ferential equations, it is valuable to have both a gen- 
eral abstract system and a specific concrete example in 
mind. Abstractly, we write the coordinates of the sys- 
tem as a single n-dimensional vector y that lives in the 
rt-dimensional phase space, and we write the equations of 
motion as a system of first-order differential equations: 



We will refer to a solution to Eq. (|2.1|l as a flow. As 




FIG. 1: The Lorenz attractor. All initial conditions except 
the origin (which is an unstable equilibrium) are attracted to 
the figure shown. 



a specific example, consider the Lorenz system of equa- 
tions: 

x = —ax + ay 

y = —xz + rx — y (2-2) 
z = xy — bz, 

where a, r, and b are constants. In the notation of 
Eq. p.lfl . we then have y = (x, y, z) and f (y) = (—ax + 
ay, —xz + rx — y,xy — bz). The Lorenz equations exhibit 
chaos for a wide variety of parameter values; in this pa- 
per, for simplicity we consider only one such set: a = 10, 
b = 8/3, and r — 28. For these parameter values, all ini- 
tial conditions except the origin asymptote to the elegant 
Lorenz attractor (Fig. P) . 



A. The deviation vector method 

The most straightforward method for calculating the 
largest Lyapunov exponent is to consider an initial 
point = y and a nearby point y^ 2) = y + Sy Q , 
and then evolve both points forward, keeping track of 
the difference 8y = y^ — y^ . If the motion is chaotic, 
then exponential separation implies that 
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so that the largest exponent is 

log [r e (t)] 



A„ 



where we write 



t 



IWII*yoll> 



(2.4) 



(2.5) 



with a subscript e that anticipates the ellipsoid axis dis- 
cussed in Sec. Ill B 31 Here || • || denotes the Euclidean 
norm (though in principle any positive-definite norm will 
do [13| ). It is convenient to display the results of this 
process graphically by plotting log [r e (i)] vs. t, which we 
refer to as a Lyapunov plot; since Eq. I|2.4|l is equivalent 
to log [r e (t)] = A max t, such plots should be approximately 
linear, with slope equal to the principal Lyapunov expo- 
nent. (In practice, to extract the slope we perform a 
least-squares fit to the simulation data, which is less sen- 
sitive to fluctuations in the value of log[r e (£)] than the 
ratio log [r e (tf)]/tf at the final time.) We refer to this 
technique as the (unrescaled) deviation vector method. 

It is important to note that, because of the problem of 
saturation, Eq. I|2.4|l does not define a true Lyapunov ex- 
ponent. In a chaotic system, any deviation <5yo, no mat- 
ter how small, will eventually saturate, i.e., it will grow 
so large that it no longer represents the local behavior 
of the dynamical system. Moreover, chaotic systems are 
bounded by definition [in order to eliminate trivial expo- 
nential separation of the form x(t) = xq e At ], so there is 
some bound B on the distance between any two trajecto- 
ries. As a result, in the infinite time limit Eq. I|2.4|l gives 



lim 



log||<fy||/ll<fyol 
t 



< lim 

t — >oo 



\ogB/\\Sy a \ 
t 



(2.6) 

In the naive unrescaled deviation vector method, the cal- 
culated exponent is always zero because of saturation. 

One solution to the saturation problem is to rescale the 
deviation once it grows too large. For example, suppose 
that we set ||5yo|| = e for some small e (say 10~ 8 ), and 
then allow the deviation to grow by at most a factor of /. 
Then, whenever ||<5y|| > / ||<5yo||, we rescale the deviation 
back to a size e and record the length Ri — ||<5y||/||<5yo| 
of the expanded vector. If we perform N such rescalings 
in the course of a calculation, the total expansion of the 
initial vector is then 



lfr/11 



N 



(2.7) 



where <5y/ is the final size of the (rescaled) separation 
vector. Applying Eq. I|2.4|) . we see that the approximate 
Lyapunov exponent satisfies 



log 



l^yoll 



N 



(2- 




FIG. 2: Comparison of the unrescaled (light) and rescaled 
(dark) deviation vector methods for calculating the principal 
Lyapunov exponent of the Lorenz system [Eq. 12.21 1. The 
slope of the rescaled line is the Lyapunov exponent (A max = 
0.905L0.003; see Sec. HID II . The initial deviation is ||<5y || = 
10 -8 , and rescaling occurs (for the rescaled method) if \\Sy\\ > 



irr 



Note the saturation of the unrescaled approach once the 



deviation has grown too large. 



The rescaled deviation vector method is not particu- 
larly robust compared to the rigorous method described 
below fSec. lIIB|) . and there are significant complications 
when applying it to constrained systems, but if imple- 
mented with care it provides a fast and accurate estimate 
for the largest Lyapunov exponent. Fig.|5]shows both the 
rescaled and unrescaled deviation vector methods applied 
to the Lorenz system [Eq. <|2.2[l ]. Note in particular the 
saturation of the unrescaled approach. We discuss the 
limitations of the rescaled method further in Sec. IIVI 



B. The Jacobian method 

Although the deviation vector method suffices for prac- 
tical calculation in many cases, in essence it amounts 
to taking a numerical derivative. For a one-dimensional 
function of one variable, we can approximate the deriva- 
tive at x — xq using 



/'(so) 



f(x +e) - f(x ) 



(2.9) 



for some e <C 1, but this prescription is notoriously in- 
accurate as a numerical calculation ^4|. Of course, it 
is better (if possible) to calculate the analytical deriva- 
tive f'(x) and evaluate it at xq. The higher-dimensional 
generalization of this is the Jacobian matrix, which de- 
scribes the local (linear) behavior of a higher-dimensional 
function. In the context of a dynamical system, this 
means that we can find the time-evolution of a small 
deviation Sy using the rigorous linearization of the equa- 
tions of motion: 



We refer to this as the (rescaled) deviation vector method. 



f(y + Sy) - f(y) = Df • Sy + 0(\\Sy\\ 2 ), 



(2.10) 
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where 



(DfV- = — 



(2.11) 



is the Jacobian matrix evaluated along the flow. For 
example, for the Lorenz system [Eq. I|2.2I) ] we have 

-a a 
Df = I r - z(t) -1 | , (2.12) 

y(t) x(t) -b 

where we write the coordinates as functions of time to 
emphasize that Eq. (|2.12l) is different at each time t. 



1. Jacobian diagnostic 

One note about Jacobian matrices is worth mention- 
ing: practical experience has shown that errors occasion- 
ally creep into the calculations leading to the Jacobian 
matrix, especially if the equations of motion are compli- 
cated. It is therefore worthwhile to note that Eq. (|2.1(J|) 
provides an invaluable diagnostic: calculate the quantity 



A = f (y + Sy) - f (y) - Df • Sy 



(2.13) 



for varying values of ||<5y||; if A does not generally scale 
as ||<5y|| 2 , then something is amiss. (The routines in |l^ 
include this important Jacobian diagnostic function.) 



2. The principal exponent 

The main value of Eq. i|2.10|l in the context of a dy- 
namical system is its combination with Eq. (|2.1|) to yield 
an equation of motion for the deviation Sy: 



f (y + Sy) 



d ( 



Sy) = f (y) 



d(Sy) 
dt l 



(2.14) 



so that (discarding terms higher than linear order) 
Eq. gives 



d(Sy) 
dt 



Df • Sy. 



(2.15) 



This equation is only approximately true for finite (that 
is, non-infinitesimal) deviations, but we can take the in- 
finitesimal limit by identifying the deviation Sy with an 
element £ in the tangent space at y. This leads to an 
exact equation for £: 



«=Df. 

dt 



(2.16) 



The initial value of £ is arbitrary, but it is convenient to 
require that ||£o|| = 1, so that the factor by which £ has 
grown at some later time t is simply ||£(t)||. 

The core of the Jacobian method for the principal Lya- 
punov exponent is to solve Eqs. i|2.1[l and i|2.16[) as a 




FIG. 3: The natural logarithm of the tangent vector length 
ri = ||£(t)|| vs. t for the Lorenz system. The slope of 
the rescaled line is the system's largest Lyapunov exponent 
(A m ax ~ 0.905). The figure and exponent are virtually iden- 
tical to the rescaled deviation method show in Fig. [5] 



coupled set of differential equations. As in Sec. Ill Al for 
chaotic systems the length of the deviation vector will 
grow exponentially, so that 



U(t)\\ 



which implies that 



An 



t 



(2.17) 



(2.18) 



For sufficiently large values of t, Eq. I|2.18|l provides an 
approximation for the largest Lyapunov exponent. It is 
essential to understand that there is no restriction on 
the length of the tangent vector £: the Jacobian method 
does not saturate. The only limitation on the size of £ 
in practice is the maximum representable floating point 
number on the computer. 



3. Ellipsoids and multiple exponents 

Although following the time-evolution of a tangent vec- 
tor £ in place of a finite deviation Sy solves the problem of 
saturation, it still only allows us to determine the prin- 
cipal exponent A max . For a system with n degrees of 
freedom, this leaves n — 1 exponents undetermined. In 
order to calculate all n exponents, we must introduce n 
tangent vectors. (We discuss the value of knowing all n 
exponents in Sec. Ill (Jl below.) Since n (linearly indepen- 
dent) vectors span an n-dimensional ellipsoid, this leads 
to a visualization of the Lyapunov exponents in terms of 
the evolution of a tangent space ellipsoid (Fig. QJ. Fig. [5] 
shows the corresponding Lyapunov plot. 

The general method is to introduce a linearly inde- 
pendent set of vectors • • • , It is con- 
venient to begin the integration with vectors that form 
the orthogonal axes of a unit ball, so that the vectors 
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of a matrix U, then Eq. (|2.19(l implies that 



FIG. 4: The Lorenz system with an evolving ellipsoid. The 
ellipsoid is calculated exactly in the tangent space (for a total 
time t — 0.4) and is superposed on the phase space for the 
purposes of visualization. There is one expanding axis (~ 
e°' 905i ) and one contracting axis (~ e -14 ' 57 *); the third axis 
has a fixed unit length (Sec. Ill D II . 
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FIG. 5: The natural logarithms of all three of the ellip- 
soid axes Ti vs. t for the Lorenz system, calculated using 
the Jacobian method (Sec. Ill Bfl . The slopes are the Lya- 
punov exponents. The three lines correspond to the expo- 
nents Ai w 0.905, A 2 w 0.0, and A 3 w -14.57 fSec. HID II . 
These values agree with the calculations in ||. 



£o 2 \ ■ ■ ■ ' £o"^} are orthonormal. Each of the se tan- 
gent vectors satisfies its own version of Eq. (|2.16|) : 



dt 



= Df ■£ 



(n) 



(2.19) 



dt 



Df U. 



(2.20) 



This equation, combined with Eq. I|2.1[) . describes the 
evolution of a unit ball into an n-dimensional ellipsoid. 

The value of the tangent space ellipsoid is this: if is 
the ith principal ellipsoid axis [and n(0) = 1], then 



i(t) 



(2.21) 



where Ai is the ith Lyapunov exponent. That is, the 
ellipsoid's axes grow (or shrink) exponentially, and if A; > 
for any i then the system is chaotic [Recall that 

we refer to the semiaxes as "axes" for brevity (Sec. HJ.] 
Turning Eq. I|2.21|l around, we can find the ith Lyapunov 
exponent by finding the average stretching (or shrinking) 
per unit time of the ith principal ellipsoid axis: 



log hit)} 



(2.22) 



If we combine the n tangent vectors to form the columns 



In practice, a more robust prescription is to record 
log [r, (t)] as a function of t and perform a least-squares 
fit to the pairs (tj,log[ri(tj)]) to find the slope A^. 

Though Eq. H2.22JI provides an estimate for the ith 
Lyapunov exponent, it requires us to find the n princi- 
pal axes of the final ellipsoid. While it is true that the 
columns of the final matrix U / necessarily span an ellip- 
soid, but they are not in general orthogonal; in particu- 
lar, the final tangent vectors do not necessarily coincide 
with the ellipsoid's principal axes. A first step in extract- 
ing these axes is to note an important theorem in linear 
algebra (see 8] for a proof): 

Theorem 1 Let A be an n x n real matrix consisting of 
n linearly independent column vectors {vi}™ =1 , and let 
{ s f}"=i be the eigenvalues and {ui}™ =1 the normalized 
eigenvectors of A T A (where A T is the transpose of A). 
Then {v^}™ =1 lie on an n-dimensional ellipsoid whose 
principal axes are {si Ui}" =1 . 

In other words, finding the principal axes of the ellip- 
soid represented by a matrix A is equivalent to finding 
the eigensystem of A T A. (We note that the ellipsoid is 
unique: any other matrix B whose columns {w"i}™ =1 lie 
on the same ellipsoid as {vi}™ =1 must necessarily give the 
same principal axes.) 

In principle, we are done: simply evolve U for a long 
time, and find the eigenvalues of U T U. In practice, this 
fails miserably; every (generic) initial vector £^ has some 
component along the direction of greatest stretching, so 
all initial tangent space vectors eventually point approx- 
imately along the longest principal axis. As a result, all 
axes but the longest one are lost due to finite floating 
point precision. 

The solution is to find new orthogonal axes as the 
system evolves. In other words, we can let the system 
evolve for some time T, stop to calculate the principal 
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axes of the evolving ellipsoid, and then continue the inte- 
gration. The method we advocate is the Gram-Schmidt 
orthogonalization procedure, which results in an orthog- 
onal set of vectors spanning the same volume as the orig- 
inal ellipsoid, and with directions that converge to the 
true ellipsoid axes. This approach, originally described 
in pj|. is a common textbook approach IIJQ, and was 
used successfully by the present author in Numeri- 
cally, the Gram-Schmidt algorithm is subject to consider- 
able roundoff error Q , and is usually considered a poor 
choice for orthogonalizing vectors, but in the context of 
dynamics its performance has proven to be astonishingly 
robust. (See Sec.|W]for further discussion.) 

We review briefly the Gram-Schmidt construction, and 
then indicate its use in calculating Lyapunov exponents. 
Given n linearly- independent vectors {u^}, the Gram- 
Schmidt procedure constructs n orthogonal vectors {v^} 
that span the same space, given by 

j=l II 3 II 

To construct the ith orthogonal vector, we take the ith 
vector from the original set and subtract off its projec- 
tions onto the previous i — 1 vectors produced by the 
procedure. The use of Gram-Schmidt in dynamics comes 
from observing that the resulting vectors approximate 
the axes of the tangent space ellipsoid. After the first 
time T, all of the vectors point mostly along the prin- 
cipal expanding direction. We may therefore pick any 
one as the first vector in the Gram-Schmidt algorithm, 
so choose £i = Ui without loss of generality. If we let 
ei denote unit vectors along the principal axes and let 
be the lengths of those axes, the dynamics of the system 
guarantees that the first vector Ui satisfies 

ui = riei + r 2 e 2 H « rjei = vi 

since ei is the direction of fastest stretching. The second 
vector V2 given by Gram-Schmidt is then 

Ui • Vi 

v 2 = Ui - - — — vi « ui - nei = r 2 e 2 , 

l|vi|| 2 

with an error of order r 2 jv\ . The procedure proceeds itcr- 
atively, with each successive Gram-Schmidt step (approx- 
imately) subtracting off the contribution due to the pre- 
vious axis direction. In principle, the system should be 
allowed to expand to a point where r 2 Cri, but (amaz- 
ingly) in practice the Gram-Schmidt procedure converges 
to accurate ellipsoid axes even when the system is orthog- 
onalized and even normalized on timescales short com- 
pared to the Lyapunov stretching timescale. As a result, 
the procedure below can be abused rather badly and still 
give accurate results (Sec. llVf) . 

4- The algorithm in detail 

We summarize here the method used to calculate all 
the Lyapunov exponents of an unconstrained dynamical 



system y = f(y) with n degrees of freedom: 

1. Construct an orthonormal matrix Uo whose 
columns (the initial tangent vectors) span a unit 
ball, and then integrate 

y = f (y) (2.24) 

and 

U = Df U (2.25) 

as a coupled set of In differential equations. We 
recommend choosing a random initial ball for 
genericity. 

2. At various times tj, replace U with the orthogonal 
axes of the ellipsoid defined by U, using the Gram- 
Schmidt orthogonalization procedure. This can be 
done either every time T, for some suitable choice 
of T, or every time the integrator takes a step. We 
have found the latter prescription to be especially 
robust in practice. 

3. If the length of any axis exceeds some very large 
value (say, near the maximum representable float- 
ing point value), normalize the ellipsoid and record 
the axis lengths 

axis at fcth rescaling) (2.26) 

at the rescaling time. Do the same if any axis is 
smaller than some very small number. 

4. Record the value of 

femax 

logrP = log + £ logif } (2.27) 

k=l 

at each time tj, where Li is the ith principal axis 
length. The second term accounts for the axis 
lengths at the fc max rescaling times. Note that if tj 
is a rescaling time itself, then log [Li(tj)] = log 1 = 
0, since by construction the ellipsoid has been nor- 
malized back to a unit ball. 

5. After reaching the final number of time steps N, 
perform a least squares fit on the pairs (tj , log r\ ) 
to find the slopes A^. Since 

log[n(t)] w Xit, (2.28) 

the slope Ai is the Lyapunov exponent correspond- 
ing to the ith principal axis. Using the Gram- 
Schmidt procedure should result in the relationship 
Ai > ... > A„. 

Most of the value of calculating A^ for i > 1 comes 
from having all n of the exponents (Sec. Ill CI below). 
Nevertheless, it is worth noting that the algorithm works 
for any value < m < n, so the method above can be 
used without alteration to find an arbitrary number of 
exponents. Fig. [S] shows the axis growth for m = 2 in 
the Lorenz system, while Fig. [S] shows the growth for 
to = n — 3. 
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FIG. 6: Close-up of Fig. showing the natural logarithms 
of the two largest ellipsoid axes vs. t for the Lorenz system, 
calculated using the Jacobian method (Sec. Ill rH . The slopes 
are the Lyapunov exponents. The plot for the larger axis 
closely matches the figures for the rescaled deviation vec- 
tor method (Fig. |5J and the single tangent vector Jacobian 
method (Fig. EH . 



C. The value of multiple exponents 

Calculating all the exponents of a system of differen- 
tial equations allows us to paint a more complete picture 
of the dynamics in several different ways. In particular, 
with all n exponents comes the ability to visualize the 
entire phase space ellipsoid (instead of just its principal 
axis), as in Fig. 0] Another important benefit of know- 
ing all the exponents is a determination of dissipative or 
conservative behavior. Conservative flows preserve phase 
space volumes, while dissipative flows contract volumes. 
Geometrically, the volume V of an ellipsoid is propor- 
tional to the product of its principal axes {r^}, so that 
the ratio of the final to the initial volume is 



(2.29) 



assuming that the initial volume is a unit ball. For dissi- 
pative systems, phase space volumes in general contract 
exponentially according to 



Yl 

V 



-M 



where A is a positive constant. Combining Eq. 
Eq. l|2~3U|> yields 



(2.30) 



and 



A 



log 



(j[rA = -^logr, = - 



the system is conservative. The special case of Hamilto- 
nian systems is of particular interest, since the equations 
of motion for many mechanical systems can be derived 
from a Hamiltonian. The Hamiltonian property strongly 
constrains the Lyapunov exponents, which must cancel 
pairwise: to each exponent +A there corresponds a sec- 
ond exponent —A |13|| . Several examples of this ±A prop- 
erty of Hamiltonian systems appear below. 

Having all the Lyapunov exponents also allows us to 
verify that there is at least one vanishing exponent, cor- 
responding to motion tangent to the flow, which must be 
the case for any chaotic system. (See Ref. @ for a proof.) 
Since we have finite numerical precision, we do not ex- 
pect to find any exponent to be identically zero, but some 
exponent should always be close to zero. A practical cri- 
terion for "close to zero" is to compute error estimates 
for the least-squares fits advocated in Sec. Ill B~4l an ex- 
ponent is "close to zero" if it is zero to within the stan- 
dard error of the fit. Applications of this method appear 
in Sec. Ill D II and Sec. Ill D 21 below. It is worth noting 
that the fitting errors are not the dominant source of 
variance in calculating Lyapunov exponents; variations 
in the initial conditions and initial deviation vectors con- 
tribute more to the uncertainty than errors in the fits. 
See Sec. IllDTl for further discussion. 

One final note deserves mention: the statement that 
A = — Xi is equivalent to a theorem due to Liou- 
ville Q , which relates the volume contraction to the trace 
of the Jacobian matrix: 



A, 



(2.31) 

where the Xi are the Lyapunov exponents. In other 
words, the phase space volume contraction constant A is 
equal to minus the sum of the Lyapunov exponents. 

If the Lyapunov exponents sum to zero, then the con- 
traction factor vanishes, and volumes are conserved — i.e, 



Yl 



exp ^ Tr Df (t) d?j , (2.32) 



where again we assume that Vq corresponds to a unit 
ball. If the trace of the Jacobian matrix happens to be 
time-independent, then this yields 
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To 



= exp [(Tr Df ) t] , (time-independent trace) 



(2.33) 

so that Eq. I|2.31|) gives A = — TrDf. In this special case, 
we can perform a consistency check by verifying that 

Xi = Tr Df . (time-independent trace) (2.34) 



D. Examples 

1. The Lorenz system 

Following the phase space ellipsoid allows us to visual- 
ize the dynamics of the Lorenz system in an unusual way. 
Fig-EJshows the Lorenz attractor together with the phase 
space ellipsoid for a short amount of time (tf = 0.4). The 
initial ball is evolved using Eq. 12.201 so it represents the 
true tangent space evolution, which is then superposed 
on the Lorenz phase space (x, y, z). It is evident that the 
initial ball is stretched in one direction and flattened in 
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another, as well as rotated. (As we shall see, the third di- 
rection is neither stretched nor squeezed, corresponding 
to the zero exponent discussed in Sec. Ill CI ) 

By recording natural logarithms of the ellipsoid axes 
as the system evolves, we can obtain numerical estimates 
for the Lyapunov exponents, as discussed in Sec. Ill B 41 
A plot of log [Tj(t)] vs. t appears in Fig. [5] for a final time 
tf = 50, with the slopes giving approximate values for 
the exponents. Using a tt = 5000 integration for greater 
accuracy yields the estimates 



Ai = 0.905 ±9 x 10~ 6 

A 2 = 1.5 x 10~ 6 ±1.7x 10 -6 

A 3 = -14.57 ±9 x 10~ 6 



(2.35) 



for the parameter values a — 10, b = 8/3, and r = 28. 
The ± values are the standard errors on the least-squares 
fit of log [ri(t)} vs. t. One of the exponents is close to zero 
(as required for a flow) in the sense of Sec. Ill Cl the error 
in the fit not small compared to the exponent. [In the 
case shown in Eq. I|2.35|l . the "error" is actually larger 
than the exponent.] The other two exponents are clearly 
nonzero, with the positive exponent indicating chaos. 

As mentioned briefly in Sec. Ill Cl the largest source of 
variance in calculating Lyapunov exponents is variations 
in the initial conditions, not errors in the least-squares 
fits used to determine the exponents. We express the 
exponents in the form 



A± 



where 



N 



N 



N ^ 



(2.36) 



(2.37) 



is the sample mean and 



\ 



N 



i N 

— T,( XU) -~ X Y 



(2.38) 



is the standard deviation. For the Lorenz system, using 
a final integration time of tf = 5000 for N — 50 ran- 
dom initial balls [all centered on the same initial value of 

(x Q ,y ,Zo)] gives 



Ai = 0.9053 ±4.1 x 10~ 4 

A 2 = -4.5 x 10~ 6 ± 7.6 x 10~ 7 

A 3 = -14.5720 ±4.1 x 10-4 



(2.39) 



The values of the error are much greater than the stan- 
dard errors associated with the least-squares fit for the 
slope for any one trial. As expected, it is evident that A 2 
is consistent with zero. 

There is a strongly expanding direction and a very 
strongly contracting direction in the Lorenz system, 
and the volume contraction constant A is large: A = 



— Ai = 13.67, so that after a time t = 5000 the vol- 
ume is an astonishingly small 6.75 x 10 -29674 . This is de- 
spite the exponential growth of the largest principal axis, 
which grows in this same time to a length 1.52 x 10 1965 ; 
the volume nevertheless contracts, since the smallest axis 
shrinks to 4.44 x 10~ 31639 in the same time. We note that 
the periodic renormalization and reorthogonalization of 
the ellipsoid axes is absolutely essential from a numer- 
ical perspective, since these axis lengths are far above 
and below the floating point (double precision) limits of 
xmax ss xmin -1 ss 10 308 on a typical IEEE-compliant 
machine [Tij . 

The Lorenz system affords an additional check on the 
numerically determined exponents: the trace of the Ja- 
cobian matrix [Eq. H2.12(l ] is time-independent, so the 
exponents should satisfy Eq. (|2.34JI : 



A, = -13.67 = TrDf = -((7+1+6) = - — » -13.67. 

(2.40) 

Eq. 12.34|l is thus well-satisfied. 

2. The forced damped pendulum 

We turn now to our second principal example of a 
chaotic dynamical system, the forced damped pendulum 
(FDP). This is a standard pendulum with damping and 
periodic forcing; written as a first-order ODE, our equa- 
tions are as follows: 



i 



p sin t 



(2.41) 



1 



Here c is the damping coefficient and p is the forcing 
amplitude, and the gravitational acceleration g and pen- 
dulum length I are set to one for simplicity. We include 
the equation t = 1 so that the system is autonomous (i.e., 
we remove the explicit time-dependence by treating time 
as a dynamical variable with unit time derivative). In 
addition to being an example with transparent physical 
relevance (in contrast to the Lorenz system), the forced 
damped pendulum, in slightly altered form, serves as a 
model constrained system in Sec. IIIII below. 

The forced damped pendulum is chaotic for many val- 
ues of c and p. For simplicity, in the present case we fix 
c = 0.1 and p = 2.5. A plot of 9 vs. t shows the system's 
erratic behavior (Fig. [7J , but a more compelling picture 
of the dynamics comes from a time-27r stroboscopic map. 
A time-T map involves taking a snapshot of the system 
every time T and then plotting uj vs. 6. Since the forcing 
term in Eq. 12.41|) is 27r-periodic, this provides a natural 
value for T in the present case. The resulting plot shows 
the characteristic folding and stretching of a fractal at- 
tractor (Fig. |8J) , which for the FDP attracts almost all 
initial conditions |8j|. 
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FIG. 7: 9 vs. t for the forced damped pendulum [Eq. 12.411 1. 




FIG. 8: lu vs. 9: the time-27r stroboscopic map for the 
forced damped pendulum. A point (u = 9, 9) is plotted ev- 
ery time 2-7T, resulting in a fractal attractor characteristic of 
dissipative chaos. 
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FIG. 9: The natural logarithms of all three of the ellipsoid 
axes n vs. t for the forced damped pendulum, calculated using 
the Jacobian method (Sec. Ill lH . The slopes are the Lyapunov 
exponents. The three lines correspond to the exponents Ai = 
0.160 ± 0.0049, A 2 = 0.0, and A 3 = -0.262 ± 0.0053 fSec.lTVl 
and Table QJ. 



The forced damped pendulum is dissipative and 
strongly chaotic. We calculate the Lyapunov exponents 
(Fig. using the Jacobian matrix: 



Df 




1 

— c p cos t 



(2.42) 











The Lyapunov exponents are (for a tf — 5 x 10 4 integra- 
tion) 



Ai 
A 2 



0.160 ±7 x 10~ 6 
8 x 10~ 8 ± 1 x 10~ 7 



(2.43) 



A, = -0.262 ± 7 x 10" 



where the error terms are the standard errors in the least- 
squares fit for the slope. (See Sec. II VI and especially Ta- 
blc[I]for the true errors due to varying initial deviations.) 
One exponent is consistent with zero (as required for a 
flow) to within the error of the fit. The dissipation con- 
stant is A = — Ai = 0.1. The trace of the Jacobian 
matrix is time-independent, so that TrDf = — c, and 
indeed A, = —0.1 = — c = TrDf as predicted by 
Eq.E33 

The zero exponent in the FDP is associated with the 
time "degree of freedom" in the Jacobian: if we delete 
the final row and column of the Jacobian matrix, only 
the positive and negative exponents remain (see, e.g., 
Fig. 1131 below). Since the time is not an actual dynam- 
ical variable, for the remainder of this paper we will 
suppress this "time piece," but it is important to note 
that the time dependence is absolutely crucial to the 
presence of chaos. According to the Poincare-Bendixon 
theorem @ , an autonomous system of differential equa- 
tions with fewer than three degrees of freedom cannot 
be chaotic. We will treat the FDP system as a time- 
dependent system with two degrees of freedom, but the 
extra equation t — 1 in the autonomous formulation is 
what creates the potential for chaos. 

An instructive case to consider is the limit c = p = 0. 
In this limit, the system is a simple pendulum, which is a 
Hamiltonian system. A simple pendulum is not chaotic, 
of course, and both its Lyapunov exponents are zero, 
but the Hamiltonian character of the system nevertheless 
shows up in the ±A property discussed above fSec. HI (J|l : 
numerically, the exponents approach zero in a symmetric 
fashion, as shown in Fig. ^3 



III. LYAPUNOV EXPONENTS IN 
CONSTRAINED FLOWS 

We come now to the raison d'etre of this paper, 
namely, the calculation of Lyapunov exponents for con- 
strained systems. For pedagogical purposes, our primary 
example is the forced damped pendulum with the po- 
sition written in Cartesian coordinates. In addition to 
this instructive example, we also discuss two constrained 
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FIG. 10: The natural logarithms of the ellipsoid axes n vs. t 
for the forced damped pendulum in the limit of zero dissi- 
pation and zero forcing (i.e., a simple pendulum). The Lya- 
punov exponents are zero, and the distance between nearby 
trajectories grows linearly (leading to logarithmic growth in 
this log plot). Nevertheless, the Hamiltonian character of the 
system is manifest in the ±A symmetry: for each exponent 
+A, there is a corresponding exponent —A. In the nonchaotic 
limiting case shown here, the Lyapunov exponents approach 
zero symmetrically. 



systems of astrophysical interest, involving the orbits of 
spinning compact objects such as neutron stars or black 
holes (see, e.g., an d EH and references therein). 

Written in terms of the Cartesian coordinates (x, y) — 
(cos 9, sin 9), the equations of motion for the FDP 
[Eq. (|2.41|) ] become (upon suppressing the time piece) 

x = —Luy 

y = lux (3-1) 
to = —cuj — y + psmt 

For a pendulum with unit radius, the Cartesian coordi- 
nates of the pendulum satisfy the constraint 

x 2 +y 2 = l. (3.2) 

Although it is certainly possible to use (x,y) in the equa- 
tions of motion, along with (x,y), this is an unnecessary 
complication; in order to keep the equations as simple 
as possible, we retain the variable to in the equations of 
motion. 

Developing the techniques for solving constrained sys- 
tems using this toy example has several advantages. The 
equations of motion and the constraint are extremely 
simple, which makes it easy to see the differences be- 
tween the constrained and unconstrained cases. In addi- 
tion, the constraint is easy to visualize, and yet it cap- 
tures the key properties of much more complicated con- 
straints. Finally, since we have already solved the same 
problem in unconstrained form, it is easy to verify that 
the techniques of this section reproduce the results from 
Sec. l!TD~2l 



A. Constraint complications 

To see how constraints complicate the calculation of 
Lyapunov exponents, consider an implementation of the 
deviation vector approach (Sec. Ill A|) . In the uncon- 
strained forced damped pendulum, given an initial condi- 
tion, we would construct a deviated trajectory separated 
by a small angle 89 (and a small velocity Slu). In the 
constrained version, a naive implementation would use a 
deviated trajectory with spatial coordinates x + Sx and 
y + Sy, where Sy = (Sx, 8y) is a small but otherwise 
arbitrary deviation vector. But the deviations are not 
independent; the deviated initial condition must satisfy 
the constraint: 

(x + Sx) 2 + (y + Sy) 2 = 1. (3.3) 

To lowest order in Sx, we must have Sy = — {x/y) Sx. 

We can now consider a more general case. Suppose 
there are k constraints, which we write as a ^-dimensional 
vector equation C(y) = 0. (In our example, C has only 
one component: with y = (x,y,u>), we have Ci(y) = 
x 2 +y 2 ~l = 0.) Then if a point y satisfies the constraints, 
the deviated trajectory must satisfy them as well: 

C(y + (5y) = 0. (3.4) 

We will refer such a as a constraint-satisfying devia- 
tion. 

Let us outline one possible method for constructing 
such a constraint-satisfying deviation. Let n be the 
number of phase space coordinates (n — 3 for the con- 
strained forced damped pendulum model). Consider 
an n-dimensional vector yo that has d nonzero entries, 
where d represents the true number of degrees of freedom 
(d = 2 for the constrained FDP). Assume that we have 
some method for constructing from yo an n-dimensional 
initial condition yo that satisfies the constraints. For ex- 
ample, we could specify the initial values of x and uj, and 
then derive an initial value of y using y = y/1 — x 2 (or 
y = — yi — x 2 ; more on this later). Now consider an n- 
dimensional vector y' = yo + Syo, which adds arbitrary 
deviations to d degrees of freedom. We can then use the 
same method as above to find y from y , and then set 

<5yo = y' - yo (3.5) 

to arrive at a constraint-satisfying deviation. 



B. Constrained deviation vectors 

Having determined <5yo by Eq. (|3.5|) (or by some other 
method), we can immediately apply the unrescaled devi- 
ation vector approach: simply track y' and y as the two 
trajectories evolve, and monitor the length of Sy = y' — y. 
Since the equations of motion preserve the constraint, 
the resulting Sy is always constraint-satisfying. The only 
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FIG. 11: Comparison of the unrescaled (light) and rescaled 
(dark) constrained deviation vector methods for calculating 
the principal Lyapunov exponent of the constrained forced 
damped pendulum (Sec. 1111 B 11 . The slope of the rescaled 
line is the Lyapunov exponent, Ai = 0.161±0.0046 (Scc. lIVl l. 
The initial deviation is ||5yo|| = 10 -6 , and rescaling occurs 
(for the rescaled method) if ||5y|| > 10~ 2 , which happens 4 
times in this figure. As in Fig. [5] the unrescaled approach 
saturates once the deviation has grown too large. 



FIG. 12: The natural logarithm of the tangent vector length 
ri = ||£(t)||r vs. t for the constrained forced damped pendu- 
lum, using a constraint-satisfying tangent vector (Sec. 1111 B 2> . 
We use the restricted norm || • || r to calculate phase space dis- 
tances (see text). Compare to Fig. [5] (unconstrained Jacobian 
method) and Fig. 1111 (constrained deviation vector method) . 



2. A Jacobian method for the largest exponent 



subtlety is using a restricted norm to eliminate the ex- 
tra degrees of freedom; for example, the restricted FDP 
norm is 

||<5y|| r = VSx 2 + 5uj 2 (3.6) 

if we choose to eliminate the y degree of freedom. Since 
Sy « —{x/y)8x, using the full Euclidean distance would 
add the term Sy 2 — (x 2 /y 2 ) Sx 2 to the expression under 
the square root, leading to an overestimate for the prin- 
cipal exponent. The restricted norm avoids this problem 
by considering only true degrees of freedom. 



1. Rescaling for constrained systems 

In contrast to the simplicity of the unrescaled method, 
the rescaled deviation vector method requires great care, 
since a carelessly rescaled deviation is not constraint- 
satisfying: C(y + Sy/r) ^ for a rescaling factor r ^ 1. 
In this case, it is necessary to extract 5y from 5y and then 
rescale it back to its initial size ||<Syo|| using the restricted 
norm. By reapplying the procedure leading to Eq. 13. 51) . 
we then find a new (rescaled) constraint-satisfying Sy 
that satisfies ||<5y||r = ||<5yo||- In this case, it is essential 
that the new deviation vector have the same constraint 
branches as the old one. For example, suppose that in 
the FDP case the value of y is negative before the rescal- 
ing. When calculating a new y' to arrive at the rescaled 
deviation <5y, it is then essential to choose the negative 
branch in the equation y' = ±yl — x' 2 . The result of 
implementing this constrained deviation vector method 
to the forced damped pendulum appears in Fig. 1111 



The method outlined above for unrescaled deviation 
vectors leads to a remarkably simple implementation of 
the single tangent vector Jacobian method. Given a 
constraint-satisfying deviation <5yo, set 

£o= Wll^ollr, (3.7) 

where || • || r is a restricted norm on the d true degrees of 
freedom. We refer to such a £ as a constraint- satisfying 
tangent vector. Since the equations of motion preserve 
the constraints, we can evolve this tangent vector using 
Eq. H2.16fl . The Jacobian method does not saturate, so 
we need only rescale if ||£|| r approaches the floating point 
limit of the computer. We can then use a procedure based 
on the rescaled deviation method to find a new (rescaled) 
constraint-satisfying tangent vector, but this is typically 
unnecessary since by the time the floating point limit 
has been reached we already have a good estimate of the 
principal Lyapunov exponent. The resulting Lyapunov 
plot for the constrained FDP appears in Fig. IT21 



3. Ellipsoid constraint complications 

We now have three methods at our disposal for calcu- 
lating the largest Lyapunov exponent, but for d degrees 
of freedom there are d exponents. What of these other 
exponents? Here we find an essential difficulty in im- 
plementing the ellipsoid method described in Sec. Ill B 31 
The core problem is this: the tangent vectors must be 
orthogonalizcd in order to extract all d principal ellip- 
soid axes, but at the same time each tangent vector must 
be constraint-satisfying. Simply put, it is impossible in 
general to satisfy the requirements of orthogonality and 
constraint satisfaction simultaneously. 
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We present here two different solutions to this problem, 
which we will refer to as the restricted Jacobian method 
and the constrained ellipsoid method. 



C. Restricted Jacobian method 

The most natural response to a system with more co- 
ordinates n than degrees of freedom d is to eliminate the 
spurious degrees of freedom using the constraints. Un- 
fortunately, this procedure is often difficult in practice: 
solving the constraint equations may involve polynomial 
or transcendental equations that have no simple closed 
form. Even for the simple case of the FDP, the sign 
ambiguity in y = — x 2 makes a simple variable sub- 
stitution impractical. Fortunately, such substitutions are 
unnecessary: since the equations of motion preserve the 
constraints, there is no need in general to eliminate n — d 
coordinates. In fact, constraints can be a virtue, since 
they can be used to check the accuracy of the integra- 
tion. 

The same cannot be said of the Jacobian matrix. As 
argued above, the extra degrees of freedom lead to funda- 
mental difficulties in applying the Jacobian method for 
finding Lyapunov exponents; constraints, far from be- 
ing a virtue, are a considerable complication. In con- 
trast to the equations of motion, though, it is relatively 
straightforward to eliminate the spurious degrees of free- 
dom. The trick is to write a restricted d x d Jacobian 
matrix, with entries only for d coordinates. 

An example should make this clear. For the FDP sys- 
tem in constrained form, we wish to eliminate one degree 
of freedom in the Jacobian matrix, and we can choose to 
eliminate either x or y. Choosing the latter, the Jacobian 
becomes 



Df 



r):r 



i).'c 



dx duj 



(3.8) 



where we have suppressed the derivatives with respect to 
the 



; 'time degree of freedom" (as discussed in Sec. lIID 21) . 
The term to focus on here is dx/dx, which seems to be 
zero a priori since x — —toy, but this is only true if we 
treat x and y as independent. Since we are eliminating 
the y degree of freedom, we cannot treat them as inde- 
pendent; y has a nonzero derivative with respect to x, so 
that 



dx 
dx 



dy_ 

dx 



(3.9) 



If we find dy/dx using y = ±vl — x 2 , we have exactly 
the same sign ambiguity problem that we had in trying 
to eliminate the y degree of freedom in the equations 
of motion. The difference here is the we need only the 
derivative of y, not an explicit solution for y in terms of x, 
and this we can achieve by differentiating the constraint: 



d_ 

dx 



(x 2 + y 2 ) = 2x + 2y 



dy_ 

dx 



dy_ 

dx 



(3.10) 
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FIG. 13: The natural logarithms of both ellipsoid axes for 
the constrained forced damped pendulum, calculated using 
the restricted Jacobian method (Sec. UllTl . The slopes are 
the Lyapunov exponents. The results agree well with the 
unconstrained case (Fig. [!|] and Tabic QJ. 



If we integrate the equations of motion using the vari- 
ables (x,y,u)), then we have the value of y at any par- 
ticular time, and we never need deal with the sign ambi- 
guity. Using the same trick to calculate duj/dx, we can 
write the restricted Jacobian as 



Df 



y 

X 

y 



y 



(3.11) 



We now proceed exactly as in the unconstrained Jaco- 
bian method, using the restricted Jacobian to calculate 
the evolution of the initial tangent space ball. Since we 
deal only with a number of coordinates equal to the true 
number of degrees of freedom, the constraints are not 
a consideration, and we can reorthogonalize exactly as 
before. 

The general case is virtually the same. For n coordi- 
nates and d degrees of freedom, there must be m = n — d 
constraint equations of the form 



C fc (y) = 



(3.12) 



for k — 1 ... to. We must choose which d coordinates 
to keep in the Jacobian matrix, eliminating to coordi- 
nates in the process. By differentiating the constraints, 
we arrive at to linear equations for the derivatives of the 
to eliminated coordinates in terms of the n variables: 



dCk 
9Vj 



= 0, 



(3.13) 



where j ranges over the indices of the eliminated coor- 
dinates (j = 2, corresponding to y, for the FDP). Since 
these are linear equations, they are both easy to solve 
and do not suffer from any sign or branch ambiguities. 
The dx d Jacobian matrices that result allow the calcu- 
lation of Lyapunov exponents with all the robustness of 
the Jacobian method for unconstrained systems. 
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FIG. 14: The orbit of a spinning relativistic binary, calculated 
using the post-Newtonian equations of motion. The equations 
model two spinning bodies, but we use an effective one-body 
approach to reduce the dynamics to the motion of one body. 
Distances are measured in terms of GM/c 2 , where M — mi + 
m2 is the total mass of the system. For a pair of black holes, 
each with 10 times the mass of the Sun, the length unit is 
GM/c 2 = 20 GMq/c 2 = 30 km. 



We considered the constrained forced damped pendu- 
lum for purposes of illustration, but it is admittedly ar- 
tificial. A more realistic example is shown in Fig. 1141 
which illustrates the dynamics of two spinning black holes 
with comparable masses. (Such systems are of consider- 
able interest for ground-based gravitational wave detec- 
tors such as the LIGO project.) The equations of motion 
come from the Post-Newtonian (PN) expansion of full 
general relativity — essentially, a series expansion in the 
dimensionless velocity v/c, where the first term is ordi- 
nary Newtonian gravity and the higher-order terms are 
post-Newtonian corrections (see, e.g., 0,^3,01)- The 
constraint comes from the spins of the black holes: it is 
most natural to think of the spin as having two degrees 
of freedom (a fixed magnitude with two variable angles 
specifying the location on a sphere), but the equations 
of motion use all three components of each hole's spin. 
We apply the methods described above to eliminate one 
of the spin degrees of freedom for each black hole, using 
the constraints 



Sli + Sl t + Si i = S, a = const., i e {1, 2}. (3.14) 

Using the effective one-body approach a priori the 
system has 12 degrees of freedom: three each for rela- 
tive position x, momentum p, and the spins Si and S2. 
Eliminating two spin components leaves 10 true degrees 
of freedom. As a result, the system has 10 Lyapunov ex- 
ponents, as shown in Fig. 1151 note in particular the ±A 
symmetry characteristic of Hamiltonian systems. 
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FIG. 15: The natural logarithms of the ellipsoid axes n vs. t 
for the system shown in Fig. 1141 Time is measured in units 
of GM/c 3 , where M = mi + 7712 is the total mass of the 
system. For two 10 solar-mass black holes, the time unit is 
GM/c 3 = 20 GM Q /c 3 = 10" 4 s. The spin magnitudes are 
fixed, so that each spin vector represents only two true de- 
grees of freedom. We deal with this constraint by using the 
restricted Jacobian method (Sec. llll"Cl . Two nonzero expo- 
nents are clearly visible, but all the others are consistent with 
zero. Note the ±A symmetry characteristic of Hamiltonian 
systems. 
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FIG. 16: The natural logarithms of both ellipsoid axes r, vs. t 
for the constrained forced damped pendulum, calculated us- 
ing the constrained ellipsoid method (Sec. 111113)1 . The slopes 
are the Lyapunov exponents. The results agree well with the 
unconstrained case (Fig. El and Table GJ. 



D. Constrained ellipsoid method 

The restricted Jacobian method relies on eliminat- 
ing spurious degrees of freedom from the Jacobian ma- 
trix, but such a prescription relies on making a choice — 
namely, which coordinates to eliminate. Each choice re- 
sults in a different Jacobian matrix. Since calculating the 
Jacobian matrix even once can be a formidable task for 
sufficiently complicated systems, it is valuable to have 
a method that uses the full Jacobian — treating all co- 
ordinates as independent — which can be calculated once 
and then never touched again. This requirement leads 
to the constrained ellipsoid method, which uses the full 
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Jacobian matrix to evolve constraint- satisfying tangent 
vectors, collectively referred to as a "constrained ellip- 
soid." When recording ellipsoid axis growth, we extract 
from each vector a number of components equal to the 
true number of degrees of freedom, resulting in vectors 
that can be orthogonalized and (if necessary) normalized 
just as in the unconstrained case. 

A detailed description of the constrained ellipsoid al- 
gorithm appears below, but we first present an impor- 
tant prerequisite: calculating constraint-satisfying tan- 
gent vectors. Let a tilde denote a vector with dimen- 
sion d equal to the true number of degrees of freedom 
(as in Sec. IIII At . We construct a full tangent vector £ 
(with n components) from a d-dimensional vector £ at a 
point y on the flow as follows: 

1. Let y' = y + e£ for a suitable choice of e. 

2. Fill in the missing components of y' using the con- 
straints to form y' as in Sec. IIII Al 



As before, we use the constrained FDP model for pur- 
poses of illustration. Treating each coordinate as inde- 
pendent yields [upon differentiation of Eq. (|3.1[) ] : 



3. Infer the full tangent vector £ using 



(3.15) 



Setting the initial conditions is now simple: form a ran- 
dom d x d matrix, orthonormalize it, and then infer the 
full d x n matrix using the method above on each column. 
The construction of constraint-satisfying tangent vectors 
described above is also necessary in the reorthogonaliza- 
tion steps of the constrained ellipsoid method. 

The full method is an adaptation of the Jacobian 
method from Sec. Ill B 41 

1. Construct a random d x d matrix and orthonor- 
malize it to form a unit ball. Use the constraints 
to infer the full d x n matrix U. 

2. Evolve the system forward using the equations of 
motion and the evolution equation for U, 



U = Df U. 



(3.16) 



3. At each time T, extract the relevant eight compo- 
nents from each tangent vector to form a d x d ellip- 
soid, orthonormalize it, and then fill in the missing 
components using the constraints, yielding again a 
d x n matrix. The restricted norms of the d tangent 
vectors contribute to the running sum for the logs 
of the ellipsoid axes [Eq. I|2.27|l ] . 

It is important to note that, unlike the other Jaco- 
bian methods, rescaling every time time T (or some 
similar method) is required for the inference equation 
[Eq. (|3.15|) ]. since the product of e and the components 
of £ must be small for the inference to work correctly. 
The method only works if the system is renormalized 
regularly, so the value of T should be chosen to be small 
enough that no principal ellipsoid axis grows too large. 
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(3.17) 



The coordinates are not independent, of course, but this 
Jacobian matrix satisfies Eq. (|2.10() as long as the de- 
viation is constraint-satisfying. For example, using the 
full deviation vector <5y = (5x, Sy, Slu) with Eq. (|3.17(l 
gives the same result as the restricted deviation vector 
<5y = (Sx, Suj) with Eq. I|3.11|) . as long as Sy — -~{x/y) Sx. 
As a result, the Lyapunov exponents calculated with the 
constrained ellipsoid method (Fig. [TB|l agree closely with 
the restricted Jacobian method (and with the original 
unconstrained results [Fig. ©])■ 

As a final example of the constrained ellipsoid method, 
consider Fig. El which shows a solution to equations 
that model a relativistic spinning test particle (e.g., a 
black hole or neutron star) orbiting a supermassive ro- 
tating black hole. (The case illustrated is a limiting case 
of the equations, which is mathematically valid but not 
physically realizable; see [l(|.) These equations (usually 
called the Papapetrou equations) are highly constrained, 
so a naive calculation of the Lyapunov exponents is not 
correct. It was the complicated nature of the Jacobian 
matrix for this system that originally motivated the de- 
velopment of the methods in this section ^j- A Lya- 
punov plot corresponding to the orbit in Fig. 1171 is shown 
in Fig. El Note especially the ±A symmetry, a result of 
the Hamiltonian nature of the equations of motion. 



IV. COMPARING THE METHODS 

A summary plot of all the methods discussed in this pa- 
per, applied to the forced damped pendulum, appears in 
Fig. El K is evident that all the methods agree closely. A 
more quantitative comparison appears in Tabled which 
gives error estimates based on integrations using fixed ini- 
tial conditions and random initial deviations. This table 
was produced by using an initial point produced from the 
final values of a previous long integration, which avoids 
any transient effects due to starting at a point not on the 
attractor. The estimates for the exponents use a final 
time of tf — 10 4 , with 100 randomly chosen values for 
the deviation vector or initial ball. All the methods agree 
on the mean exponents to within one standard deviation 
of the mean. (Recall that we omit the zero exponent 
associated with the time "degree of freedom.") 



A. Speed 

The various methods for calculating the exponents dif- 
fer significantly in their execution time, as shown in Ta- 
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FIG. 17: The orbit of a small spinning compact object (such as a solar-mass black hole) in the spacetime of a rotating 
supermassive black hole, (a) The orbit embedded in spherical polar coordinates; (b) the orbit's projection onto the x-y plane. 
The lengths are expressed in terms of GM/c 2 , where M is the mass of the central black hole. For a maximally spinning black 
hole, the horizon radius is m = GM/c 2 . For the supermassive black hole at the center of the Milky Way, M = 3 x 10 6 Mq [l9||. 
which corresponds to a length unit of GM/c 2 = 4.4 x 10 9 m. The system shown here is chaotic fFig. 1181 . although this orbit 
represents a limiting case of the equations that is not physically realizable fTojl . 



TABLE I: Comparison of different Lyapunov exponent methods applied to the forced damped pendulum. We consider both 
the unconstrained [Eq. 12.411 1 and constrained [Eq. 13. H i formulations. The integrations have a final time tj = 10 4 , and for 
each method we consider 100 random initial deviations. We calculate the positive exponent (Ai) and, if possible, the negative 
exponent (A3) as well. (We omit the zero exponent (A2) for brevity.) The error estimates are the standard deviations in the 
mean, a/yN. The deviation vector methods are all rescaled. The constrained ellipsoid method rescales and reorthogonalizes 
every time T = 1, and uses a value of e — 10 -6 for the tangent-vector inference [Eq. (I3.15H 1. The error goal is a fractional error 
of 10 -10 per step. 



Method 


Ai 


A3 


unconstrained deviation vector 
unconstrained Jacobian 
constrained deviation vector 
constrained Jac. (1 tangent vector) 
restricted Jacobian 
constrained ellipsoid 


0.1610 ±0.00050 
0.1608 ± 0.00050 
0.1608 ± 0.00051 
0.1605 ±0.00048 
0.1607 ±0.00048 
0.1605 ±0.00050 


-0.2618 ±0.00053 

-0.2614 ±0.00055 
-0.2617 ±0.00051 



blc [n] Generally speaking, the deviation methods are 
faster than their Jacobian method counterparts, which is 
no surprise — the deviation vector methods involve fewer 
differential equations. More surprising is the perfor- 
mance penalty for the restricted Jacobian method. This 
is the result of a significantly smaller typical step-size 
in the adaptive integrator needed to achieve a particular 
error tolerance. The restricted Jacobian may result in 
a system of equations that is more difficult to integrate 
because of the elimination of simple degrees of freedom 
with the potentially complicated solutions to the con- 
straint derivative equations dCk/dyi — [Eq. Q3.13)) ]. 
On the other hand, the performance penalty of the re- 
stricted Jacobian method is probably worth the gain in 
robustness, as discussed below. Moreover, for other sys- 
tems (e.g., the system shown in Figs. I14land ll5fl . the re- 
stricted Jacobian method is comparable in speed to the 



other Jacobian methods. 



B. Robustness 

Numerical methods are more useful if they are rela- 
tively insensitive to small changes in implementation de- 
tails, and the Jacobian methods win in this category. 
When reorthogonalization occurs every time step, with- 
out rescaling, the plain Jacobian method is virtually bul- 
letproof. The rescaling in this case can even occur only 
when the tangent vector norms reach very large or small 
values, say ||£|| w 10 ±10 °. This robustness also applies 
to the restricted Jacobian method, which is considerably 
less finicky than any other method for constrained sys- 
tems, and we recommend its implementation if practical. 

Jacobian methods that rescale and reorthogonalize ev- 
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FIG. 18: The natural logarithms of the ellipsoid axes for 
the system shown in Fig. 1171 vs. relativistic proper time r, 
in units of GM/c s , where M is the black hole's mass. For 
the supermassive black hole at the center of the Milky Way, 
M = 3 x 10 6 Mq jTgfl . which corresponds to a time unit of 
GM/c 3 — 15s. The largest Lyapunov exponent is A ma x ~ 5x 
10~ 3 (GM/c 3 ) -1 , which corresponds to an e-folding timescale 
of t a = 1/A = 2 x 10 2 GM/c*. For M = 3 X 10 6 M Q , this 
means that nearby trajectories diverge by a factor of e in 
the local (Lorentz) frame of an observer on this orbit in a 
time r = 3000 s = 50min. We find nonzero exponents in 
this system only for physically unrealistic values of the small 
body's spin (|11|). 
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FIG. 19: Natural logarithms of the ellipsoid axes vs. t for the 
unconstrained deviation vector method (dashed), the uncon- 
strained Jacobian method from Fig. |§] (thick) and all the con- 
strained methods. The constrained methods include the fol- 
lowing: rescaled deviation vector (black), Jacobian with single 
constraint-satisfying tangent vector (red) , restricted Jacobian 
(orange), and constrained ellipsoid (dashed blue). (The colors 
appear as shades of gray in print versions of this paper.) All 
the constrained methods start with exactly the same initial 
conditions. 



TABLE II: Timing comparison for different Lyapunov expo- 
nent methods applied to the forced damped pendulum. The 
times (on a 2 GHz Pentium 4) for a final time of tf — 10 4 
are in seconds: ti for the positive exponent Ai and tis for 
the negative exponent A2; we omit the zero exponent (A2) 
for brevity. (We write 1—3 to emphasize that calculating A3 
also calculates Ai as a side-effect.) We consider both the un- 
constrained [Eq. 12.411 1 and constrained [Eq. 13. H i formula- 
tions. The integrations use a C++ Bulirsch-Stoer integrator 
adapted from (l4| . The deviation vector methods are rescaled, 
and the constrained ellipsoid method rescales and reorthogo- 
nalizes every time T = 1. The error goal is a fractional error 
of 10~ 10 per step. The relatively small difference between 
deviation vector and Jacobian methods is the result of the 
small number of degrees of freedom; for larger systems (with 
larger Jacobians) the difference can become large |Tljl . We 
note that the restricted Jacobian method is unusually slow 
for the forced damped pendulum, but this is not generally 
the case. 



Method 


ti 


h-3 


unconstrained deviation vector 


2.57 




unconstrained Jacobian 


3.65 


5.16 


constrained deviation vector 


3.51 




constrained Jacobian (1 tangent vector) 


4.05 




restricted Jacobian 


35.3 


45.0 


constrained ellipsoid 


4.30 


5.88 




t 



FIG. 20: The natural logarithms of the two larger ellipsoid 
axes for the Lorenz system using the Gram-Schmidt algo- 
rithm, with the axes rescaled every T = 10 -3 . The largest 
and smallest directions differ by less than 2% when rescaling 
this frequently, but the axes nevertheless converge rapidly to 
the correct directions (as determined by the Jacobian method, 
Fig-inj. Numerical investigations confirm that this robustness 
persists at least down to T — 10 -5 . 



ery time T are less robust, since a priori we have no 
knowledge of appropriate values for T. Experimentation 
in this case is required to find good values of T; for the 
Lorenz system, T — 1 works well, but T = 5 leads to 
inaccurate estimates for the negative exponent, as seen 
in Fig. |^ It is better to err in the direction of small 
times, since the Gram-Schmidt procedure is quite robust: 



even when rescaling occurs on very short timescales — so 
that the longest axis has almost no chance to outgrow 
the other principal axes — the Gram-Schmidt method still 
converges to the correct exponents (Fig. 120(1 . Using the 
Gram-Schmidt algorithm to find the principal axes ben- 
efits from a strong feedback mechanism, insuring rapid 
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FIG. 21: The natural logarithms of the ellipsoid axes for 
the Lorenz system, with reorthogonalization/rescaling every 
time T = 1 (dark dots) and T = 5 (light lines). The two 
larger exponents agree exactly, but the negative exponent is 
incorrect due to roundoff error, since the smallest axis shrinks 
from unity to a size of e - 5x14 - 57 ~ 2 x 1CT 32 in a time T = 5. 



convergence to the correct axes. Using a very small value 
for T greatly increases the execution time, of course. A 
useful prescription in practice is to do a short integration 
with T chosen to be small compared to any characteris- 
tic timescales in the problem, in order to obtain a first 
estimate for the exponents. We may then choose T to 
be as large as we like, consistent with the avoidance of 
unacceptable roundoff error. 

The constrained ellipsoid method is dependent on fre- 
quent rescaling to keep the size of the tangent vec- 
tors small, since the inference scheme represented by 
Eq. (|3.15|) fails for large vector norms. As a result, this 
method suffers from the complexity of all time T meth- 
ods, i.e., it requires care in choosing an appropriate value 
of T. In addition, the value of e in Eq. (|3.15|) must be 
chosen carefully to achieve accurate tangent-vector infer- 
ences: the method relies on small values of e for accu- 
racy, but values that are too small suffer from roundoff 
errors. It is advisable to calibrate the value of e so that 
the largest Lyapunov exponent agrees with the result of a 
second method (such as the single tangent-vector method 
or the deviation vector method), as discussed in [Io| . 
Such a calibration was required to produce the values 
in TableQ] the largest exponent calculated using the con- 
strained ellipsoid method differs from the other methods 
by several standard deviations when using e = 1CP 5 for 
the inference, but agrees well when using e = 10~ 6 . 

Finally, the deviation vector methods are all very fast, 
but they are sensitive to the size eo of the initial deviation 
vector. The rescaled methods are particularly inaccurate 
if the value of eo is too small, which leads to roundoff error 
in the initial size of the deviation vector and can give 
inaccurate results, as shown in Fig. [23 These methods 
should be used with care, and should always be double- 
checked with a Jacobian method if possible. 



FIG. 22: The natural logarithms of the largest ellipsoid axis 
for the constrained forced damped pendulum, calculated us- 
ing the rescaled deviation vector method for varying sizes of 
the initial deviation. We vary the size of the initial deviation 
vector from eo = 10 -4 (bottom) to eo = 10~ 13 (top). Values 
of eo between 10 -4 and 10 -8 agree closely, but smaller values 
lead to erroneously high values for the Lyapunov exponent. 
It is important to calibrate the deviation vector method using 
the Jacobian method (Sec. Ill HI if possible. 



V. SUMMARY AND CONCLUSION 

Chaotic solutions exist for an enormous variety of non- 
linear dynamical systems. Lyapunov exponents provide 
an important quantitative measure of this chaos. We 
have presented a variety of different methods for calcu- 
lating these exponents numerically, both for constrained 
and unconstrained systems. Both types of systems can be 
investigated using deviation vector methods or Jacobian 
methods. Deviation vector methods use the equations of 
motion to evolve two nearby trajectories in phase space 
to determine the time-evolution of the small deviation 
vector joining the trajectories. This family of methods 
is computationally fast, but yields only the largest expo- 
nents, and also suffers from sensitivity to the size of the 
initial deviations. The Jacobian methods share the use of 
the Jacobian matrix of the system as a rigorous measure 
of the local phase-space behavior. They are computa- 
tionally robust in general, and can be used to determine 
multiple exponents, but this comes at the cost of execu- 
tion speed. 

Calculating Lyapunov exponents for constrained sys- 
tems presents a variety of complications, all revolving 
around the notion of constraint-satisfying deviations: 
"nearby" trajectories must be chosen carefully to insure 
that they satisfy the constraints. We have presented 
several methods for dealing with these complications, 
including a deviation vector method and two Jacobian 
methods: the restricted Jacobian method, which elimi- 
nates spurious degrees of freedom in the Jacobian by dif- 
ferentiating the constraints; and the constrained ellipsoid 
method, which uses the full Jacobian matrix to evolve 
constraint-satisfying tangent vectors. These methods al- 
low the determination of all d Lyapunov exponents for 
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systems with d degrees of freedom. 
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APPENDIX A: ELLIPSOID AXES AND THE 
SINGULAR VALUE DECOMPOSITION 

In this appendix, we discuss an alternative method 
for calculating the ellipsoid axes used in the Jacobian 
method, namely, calculating the ellipsoid axes exactly. 
The method described seems superior on paper to the 
Gram-Schmidt technique described in Sec. Ill Bl but suf- 
fers from subtle complications that make it fragile in 
practice. Nevertheless, within a narrow range of va- 
lidity (specified below), calculating exact ellipsoid axes 
provides valuable corroboration of the principal Jacobian 
method discussed above. 

Recall Theorem ^ from Sec. Ill B 31 which relates the 
eigensystem of the matrix A T A to the ellipsoid spanned 
by the columns of A. In order to find the axes of an 
evolving ellipsoid, we could apply Theorem^directly, but 
there is a mathematically equivalent prescription that 
is numerically virtually bulletproof, namely, the famous 
singular value decomposition: 

Theorem 2 Let A be a nonsingular n x n matrix. Then 
there exist orthonormal n x n matrices U and V , and a 
diagonal matrix S, such that 

A = USV T . (Al) 

This is the singular value decomposition (SVD) of A, 
and the values Si in S — diag(si, . . . , s n ) are the singular 
values . 

Since V is an orthogonal matrix, we have V T = V -1 , so 
that Eq. IjAlfl is equivalent to AV = US. Geometrically, 
this means that the image of the unit ball V is equal to an 
ellipsoid whose ith principal axis is given by Si times the 
ith column of U. V in this context is a special ball, but 
the image of any unit ball is the same unique ellipsoid. 
This leads to the following theorem: 

Theorem 3 Let A be a nonsingular n x n matrix, and 
let U and S be the matrices resulting from the singular 
value decomposition of A [Eq. \A Then the columns 
of A span an ellipsoid whose ith principal axis is Sj Uj ; 
where S = diag(si, . . . , s n ) and {u^}" =1 are the columns 
ofU. 

We thus see that the singular value decomposition is 
equivalent to finding the eigensystem of A T A. (See Ap- 
pendix A in |8| for proofs of these theorems.) 
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FIG. 23: The natural logarithms of the two larger ellipsoid 
axes for the Lorenz system using the singular value decompo- 
sition. The axes are rescaled every T — 0.5 to exaggerate the 
deviations from the correct results, but any rescaling causes 
the SVD method to fail (see text). Compare to unrescaled 
SVD (Fig. 12 II and the Gram-Schmidt method with frequent 
rescaling (Fig. I2UB . 
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FIG. 24: The natural logarithms of the two larger ellipsoid 
axes for the Lorenz system using the singular value decom- 
position, without rescaling. The results agree well with the 
Gram-Schmidt method (Fig. Rescaling (which is always 
necessary if we approach the floating point limits of ~10 ±308 ) 
ruins the agreement. 

Substituting the singular value decomposition for the 
Gram-Schmidt procedure leads to a replacement of 
step (2) from Sec. ITlBl 

(2') At various times tj, replace U with the orthogonal 
axes of the ellipsoid defined by U, using the singu- 
lar value decomposition. This can be done either 
every time T, for some suitable choice of T, or ev- 
ery time the integrator takes a step. It is essential 
to order the principal axes consistently. We recom- 
mend sorting the axes so that si > S2 > . . . > s„. 

Unfortunately, this prescription behaves badly when 
rescaling is necessary, as shown in Fig. [23 The underly- 
ing cause of this is a fundamental property of the singular 
value decomposition: it is only unique up to a permuta- 
tion of the ellipsoid axes. If we adopt an ordering based 
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on the axis lengths, we can refer, for example, to the 
longest axis as axis 1. During any particular time pe- 
riod, axis 1 may grow or shrink; the only requirement 
is that it be the fast-growing axis on average. Unfortu- 
nately, rescaling the axes causes this ordering method to 
fail: if axis 1 should happen to contract between rescaling 
times, then the ordering based on length leads to incor- 
rect axis labels, since axis 1 is no longer the longest axis. 
Even worse, when ordering by axis length, the length of 
the longest axis is always added to the running sum for 
the largest Lyapunov exponent, while the length of the 
smallest axis always contributes to the smallest exponent. 
This selection bias leads to systematic errors, guarantee- 
ing overestimates for the absolutes values of both the 
exponents (Fig. I2HJI. 

If the system is not rescaled, there is still some initial 
ambiguity in axis labels, but once axis 1 has grown suf- 



ficiently large it is very unlikely ever to become smaller 
than the other axes. Thus, after an initial expansion and 
contraction phase that establishes the ordering, the axis 
labels remain fixed, and the results of the (unrescaled) 
SVD method agree well with Gram-Schmidt (Fig.[^}. 

It should be possible in principle to follow the axis 
evolution by tracking the continuous deformation of the 
ellipsoid. This would mean assigning labels to the axes 
and then ensuring, e.g., that axis 1 at a later time is in- 
deed the image of the original axis 1. This method would 
require following the system over very short timescales 
to guarantee the correct tracking of axes, and even then 
is likely to be fragile and error-prone. Because of these 
complications, we recommend the simpler Gram-Schmidt 
process, which has proven to be reliable and robust in 
practice. 
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